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ABSTRACT 

A  few  years  ago,  the  class  of  Essentially  Noii-Oscillatory  Schemes  for  the  numerical  sim¬ 
ulation  of  hyperbolic  equations  and  systems  was  constructed.  Since  then,  some  extensions 
have  been  made  to  multidimensional  simulations  of  compressible  flows,  mainly  in  the  context 
of  very  regular  structured  meshes.  In  this  paper,  we  first  recall  and  improve  the  results  of  an 
earlier  paper  about  non-oscillatory  reconstruction  on  unstructured  meshes,  emphasising  the 
effective  calculation  of  the  reconstruction.  Then  we  describe  a  class  of  numerical  schemes 
on  unstructured  meshes  and  give  some  applications  for  its  third  order  version.  This  demon¬ 
strates  that  a  higher  order  of  accuracy  is  indeed  obtained,  even  on  very  irregular  meshes. 


'Research  was  supported  by  the  National  Aeronautics  and  .Space  Administration  under  NASA  C'ontract 
No.  NASI- 19480  while  the  author  was  in  residence  at  the  Institute  for  C'omputer  Applications  in  Science 
and  Engineering(IC'ASE),  NASA  Langley  Research  (Center,  Hampton,  VA  2.4681-0001. 


1  Introduction 


During  the  past  few  years,  a  growing  interest  has  emerged  in  building  high  order  accurate 
and  robust  schemes  for  compressible  flow  simulation.  One  of  the  difficulties  is  the  appearance 
of  possibly  strong  discontinuities  that  may  interact  together,  even  for  smooth  initial  data. 
One  way  to  avoid  this  difficulty  is  to  use  a  Totally  Variation  Diminishing  scheme.  Such 
a  scheme  has  the  property,  at  least  for  ID  scalar  equations,  of  not  creating  new  extrema, 
and  hence  provides  a  reasonable  treatment  of  discontinuities.  TVD  schemes  have  since  been 
successfully  and  widely  used  with  many  types  of  meshes  {see  for  example  [1]  for  a  review, 
and  among  many  others  [2]  for  simulations  on  finite  elements  type  meshes).  Nevertheless, 
one  of  their  main  weakness  is  that  the  order  of  accuracy  falls  to  first  order  in  regions  of 
discontinuity  and  at  extrema,  leading  to  excessive  numerical  dissipation. 

Various  methods  have  been  proposed  to  overcome  this  difficulty  (for  example,  mesh 
adaptation  [3,  4,  5]),  but  one  promising  approach  may  be  the  class  of  Essentially  Non- 
Oscillatory  schemes  (E.N.O.,  for  short),  introduced  by  Harten,  Osher  and  others  [6,  7,  8, 
9,  10].  The  basic  idea  of  E.N.O.  schemes  is  to  use  a  Lagrange  type  interpolation  with  an 
adapted  stencil:  when  a  discontinuity  is  detected,  the  procedure  looks  for  the  region  around 
this  discontinuity  where  the  function  is  the  smoothest.  This  reconstruction  technique  may 
be  applied  either  to  the  node  values  [9]  or  to  specially  constructed  averages  over  control 
volumes  [6,  7,  8].  In  the  latter  case,  the  approximation  is  done  in  a  conservative  way.  This 
enables  approximation  of  any  piecewise  smooth  function  to  any  desired  order  of  accuracy. 

Until  now,  very  few  attempts  have  been  made  to  adapt  these  ideas  to  multidimensional 
flows  (see  for  example  [9,  11]),  to  smoothly  varying  grids,  or  especially  to  unstructured  grids. 
For  the  latter  topic,  only  preliminary  work  exists  (see  [12,  13,  14,  15]).  In  [14]  general 
ideas  were  presented  with  a  review  of  the  existing  ENO  methods,  but  no  implementation 
was  given.  In  [15]  several  algorithms  were  presented  and  tested  on  simple  problems  (linear 
advection.  Burger’s  like  equations),  but  their  reconstruction  algorithms  appears  to  be  very 
complicated  and  costly.  There  was  no  study  of  the  numerical  stability  of  the  reconstruc¬ 
tion.  Our  experience  has  shown  that  this  point  is  fundamental.  In  [12]  we  studied  two 
reconstruction  methods  based  on  two  different  polynomial  approximations,  and  have  also 
estimated  the  behavior  of  their  leading  coefficients.  This  enabled  us  to  design  an  algorithm 
that  has  been  shown  to  give  third  and  fourth  order  approximation.  We  also  discussed  the 
choice  of  candidate  stencils:  a  few  stencils  suffice.  In  [13]  this  reconstruction  method  was 
implemented  for  compressible  flows  problems  and  tested  on  a  2D  shock  tube  problem  on  a 
triangular  mesh.  In  the  finite  volume  scheme  used,  the  control  volumes  were  the  triangles 
of  the  mesh.  If  one  wants  this  set  to  be  as  isotropic  as  possible,  the  minimum  numl)er  of 
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possible  stencils  is  much  larger  than  in  the  version  presented  in  [12];  its  construction  is  also 
less  natural.  This  makes  the  scheme  of  [13]  very  costly. 

This  paper  is  organized  as  follows.  In  section  2  we  recall  and  improve  the  results  ob¬ 
tained  in  [12].  These  results  are  valid  whatever  the  type  of  the  underlying  mesh  (structured 
or  not).  In  particular,  we  pay  much  attention  to  the  study  of  the  numerical  stability  of  the 
reconstruction  algorithm.  In  section  3  we  present  our  numerical  scheme.  In  section  4  some 
numerical  results  are  presented.  Finally,  several  general  comments  are  made  in  the  conclu¬ 
sion.  Throughout  this  paper,  the  values  used  are  the  average  values  in  the  given  control 
volumes. 

2  The  reconstruction  problem  on  unstructured  meshes 


Let  us  first  recall  basic  facts  about  ID  reconstruction.  They  show  why  a  new  method  has 
to  be  introduced  for  unstructured  meshes.  We  first  recall  how  to  interpolate  data  in  an 
essentially  non-oscillatory  Lagrange  fashion,  and  then  how  this  is  used  to  reconstruct  ID 
data. 


Essentially  non-oscillatory  interpolation.  This  relies  on  two  well  known  properties  of 
divided  differences.  Let  {yo  <  j/i  <  .  •  -  <  J/a,}  be  a  stencil,  and  let  [i/o^  •  •  •  J/a]  the  k  +  1'^* 
divided  difference  of  v,  a  piecewise  smooth  real  valued  function. 


(i)  If  V  \s  p  >  n  times  continuously  differentiable  in  the  interval  [3/0,?/^],  then 

[?/o,  •  •  ■  ,Va,  ]  V  =  — -f  0{\yk  -  J/or"*")  for  all  k  <  n  and  some  ^  G  [yo,  y^]- 


(ii)  If  p  <  71,  admits  a  jump  [n^^']  in  [j/o,  J/t],  then 


[yo,  •  •  •  Vk]  v^O 


[n^P)] 

|yA  -  yol'^-'’ 


for  all  k,  7;  +  1  <  A:  <  n. 


These  relations  show  that  the  divided  differences  remain  bounded  whatever  the  mesh  size, 
for  smooth  functions,  but  go  to  infinity  rather  c|uickly  for  nonsmooth  functions.  With  the 
help  of  these  two  remarks,  Ilarten  et.  al.  have  derived  the  following  F.N.O.  interpolation 
algorithm:  Assume  a  grid  {,yj},?yj  <  yj+\-  For  any  j,  first  consider  =  {yj). 

(i)  If  l[.yj’.yj-f-i]”l  <  |[.yj-i,yj]»’l  then  =  {yj,yj+,},  else  =  {yj_, 
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(ii)  Assume  that  5^^^  =  {j/jg,  •••,  },  a  stencil  for  +  order  reconstruction,  is  given  (j/jo 

mdyj^  are  the  extreme  points  of  the  stencil).  U\[yjo-\, Vjo,  ,yk]v\  <  Ibjo)  •  •  • 
then  else  5^''+^)  =  5^*=) U{yi*+i }  • 

Once  the  required  number  of  points  has  been  chosen,  one  can  compute  Lagrange  interpolation 
based  on  the  last  stencil  of  the  algorithm:  this  is  the  E.N.O.  interpolation  of  v. 


The  ID  conservative  reconstruction.  We  consider  a  mesh  on  IR,  {x{}i^Zi  that  may  or 
may  not  be  regular.  Around  each  point,  a:,,  we  define  a  control  volume  [a:,_i/27  3:,q.i/2]  where 
as  usual, 

Xi  "L 


Xi+l/2  — 


Let  us  consider  u,  a  piecewise  smooth  real  valued  function.  We  let  u,  denote  the  average  of 
u  in  [xi-i/2,a:i+i/2]: 


Ui  = 


tXx-yill 

I  u{t)  dt 
Jx,-\I2 

There  are  two  classical  ways  of  reconstructing  u  from  its  averages: 


^j+l/2  “  ^j-1/2 


(1) 


(i)  Reconstruction  via  primitive  functions.  One  considers  W  a  primitive  of  u,  say 

W{x)  =  /  u(<)  dt. 

J^l/2 

The  values  of  W  at  the  points  Xj+^i2  are  easily  recovered  from  the  data: 


j 

i^(2^.+i/2)  =  12(a:.+i/2  -  a:,-i/2)tli 

j=o 


One  computes  an  essentially  non-oscillatory  reconstruction  R(W,n  +  1)  of  VE,  up  to 
the  order  n  +  1,  as  explained  above.  Here,  one  sets  yj  =  Xj+if2.  The  reconstruction  of 
u,  R(u,n)  is  defined  as: 


R{u,n) 


dR{Wxn  +  l) 
dx 


Clearly  the  average  values  of  R{u,n)  over  any  [a:,_i/2,  a:,+i/2]  is  u,. 


(ii)  Reconstruction  via  deconvolution.  Here,  the  mesh  must  be  regular,  Xi+1/2  — 
Xi-\/2  =  Ax  .  One  may  see  equation  (l)  as  the  convolution  product  n  of  it  and  the 
characteristic  function  of  [—Ax/2,  Ax/2].  One  has  v(xi)  =  u,.  Then  v  is  reconstructed 
in  an  essentially  non-oscillatory  fashion  as  explained  above,  where  yj  =  Xj.  Finally, 
one  applies  a  deconvolution  operator  to  R{v,n),  as  in  [6]  for  example,  to  get  /?(it,n). 
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Atkins  k,  Casper  [11]  have  used  a  tensor  product  of  1  D-reconstructions  to  derive  their 
numerical  scheme.  This  is  possible  because  they  assume  a  regular  transformation  between  a 
Cartesian  mesh  of  [0,  1]  x  [0,  1]  and  their  computational  grid.  The  same  trick  has  been  used 
in  Shu  et.  al.  [9].  In  the  context  of  unstructured  grid,  these  tricks  cannot  work  for  at  least 
two  reasons: 

(i)  The  reconstruction  via  deconvolution  needs  very  regular  meshes  (each  control  volume 
can  be  obtained  from  any  other  by  translation), 

(ii)  The  reconstruction  via  primitive  function  methods  needs  to  know  point  values  of  a 
primitive  over  any  rectangle  on  the  data.  Hence,  these  rectangles  must  be  exactly 
covered  by  control  volumes.  As  can  be  seen  in  Figure  1,  this  is  in  general  impossible. 

These  two  remarks  show  that  a  straightforward  extension  of  one  dimensional  ideas  is  not 
sufficient. 

2.1  Preliminaries 

In  the  sequel,  the  symbol  IR,i[A',  T]  denotes  the  set  of  polynomials  P  in  the  variables  X  and 
Y  of  total  degree  less  than  or  equal  to  n: 

p(x,Y)  =  j2  Y,  “u-’cr' 

/=!  i+j=l 

The  set  IRn[A",  T]  is  a  vector  space  of  dimension  N{n)  =  i!i±iK!i±2l^  3^  basis  of  which  is  the 
set  of  monomials 

where  (xo,j/o)  's  any  point  of  IR^.  The  total  degree  of  P  does  not  depend  on  the  choice 
of  {xo.yo).  As  we  will  show  later,  this  kind  of  basis  is  not  the  best  suited  for  practical 
calculations. 

Let  M.  be  a  mesh  of  the  finite  element  type.  Associated  with  this  mesh,  we  have  a 
triangulation  T.  We  may  consider  several  kind  of  control  volumes,  for  example  the  triangles 
of  T  themselves  or  the  dual  mesh  (.see  Figure  2).  The  dual  mesh  are  constructed  as  follows: 
for  each  vertex  Mi,  the  control  volume  is  obtained  by  connecting  the  midpoints  of  the  edges 
incident  on  it  to  the  barycenters  of  the  surrounding  triangles  to  which  it  belongs.  Let  us 
denote  by  {6’,}  the  set  of  control  volumes.  We  only  reqtiire  the  followiTig  properties: 

•  For  any  i  /  j,  Ci  P\Cj  is  of  empty  interior. 
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•  Ci  is  connected, 


•  There  is  an  algebraic  dependency  of  the  C^’s  on  the  points  of  M  (This  is  true  for  the 
two  above  examples). 


•  The  boundary  of  Ci  is  a  locally  regular  curve  (This  is  also  true  for  the  two  above 
examples). 


We  consider  the  following  problem  (problem  'P  or  approximation  in  the  mean  for  short): 


Let  u  be  a  regular  enough  function  (say  in  ).  Given  N  and  n  two  integers,  a 
set  of  control  volumes  S  =  {C„}i</</v,  find  an  elemerit  P  G  IRijA",  F]  such  that 
for  I  <  I  <  N , 


<  u  >c,^ 


def 


JC. 


•l 


U  UrU./ 


area 


(Ci,) 


=  <  P>C, 


(2) 


For  this  problem  to  have  a  unique  solution,  one  must  fulfill  two  conditions; 


•  N  =  =  N{n) 


•  The  following  Vandermonde  matrix  must  be  non  singular: 


V  = 

[<  >C.,} 

i  +  j  <  n 

l<l<  N 

<  >0., 

<  y  >c,. 

...  <a:">c,, 

..  <K">c, 

... 

-  ’/v  ’/v  -yv  ‘TV  -  j 

(3) 

If  A5  =  det  V  0,  then  we  say  that  this  stencil  is  admissible.  In  that  ceise,  there  is  a  unique 
solution  to  problem  V  that  will  be  denoted  by 

A  similar  problem  was  first  considered  by  Barth  et.  al.  [16]  for  smooth  functions,  then  by 
Harten  et.  al.  [14],  Vankeirsblick  et.  al.  [15],  and  Abgrall  [12].  In  the  three  first  references 
[16,  14,  15],  the  authors  consider  overdetermined  systems  for  two  reasons.  First,  the  problem 
V  does  not  always  have  a  unique  solution.  Second  they  claim  that  the  condition  number 
of  the  overdetermined  system  is  better  than  that  of  problem  V.  In  [12],  the  same  approach 
as  here  was  adopted.  To  support  that  choice,  one  must  notice,  as  explained  in  remark  1, 
that  (3)  is  generally  not  singular.  Also,  the  condition  number  of  the  linear  system  depends 
mainly  on  the  basis  used  for  the  polynomial  expansion,  as  shown  in  section  2.4.  For  these 
two  reasons,  we  have  preferred  this  approach,  which  also  has  the  advantage  of  simplifying 
the  coding  of  the  global  scheme. 


Remarks: 


(i)  We  do  not  know  if  the  admissibility  condition  admits  a  geometric  interpretation  (except 

for  n  =  1).  We  do  not  even  know  whether  there  is  a  systematic  way  of  constructing 
admissible  stencils,  as  is  the  case  for  the  Lagrange  interpolation  [17].  Nevertheless, 
one  may  say  that,  in  general,  any  stencil  is  admissible:  one  may  consider  the  equation 
As  =  0  as  an  algebraic  surface  in  for  some  integer  k  This  surface  is  then  of 

empty  interior,  from  a  topological  point  of  view,  so  that  if  S  is  not  admissible,  one  only 
has  to  change  slightly  the  elements  of  S  for  it  to  become  admissible.  Nevertheless,  the 
condition  number  of  the  linear  system  may  be  very  bad.  We  will  discuss  that  point  in 
section  2.4. 

(ii)  This  admissibility  condition  is  independent  of  the  basis  chosen  for  expanding  the  poly¬ 
nomial  P. 

2.2  Some  general  results  about  approximation  in  the  mean 

In  this  section,  we  give  two  results  on  reconstruction  in  the  mean  for  a  function  u,  which 
may  or  may  not  be  smooth.  These  results  generalize  well  known  properties  of  the  Lagrange 
interpolation  of  ID  real  valued  functions,  that  have  been  used  as  a  corner  stone  by  Harten 
and  his  coauthors  to  design  an  essentially  non-oscillatory  reconstruction.  We  have  to  empha¬ 
size  that  the  reconstruction  is  the  mean  in  7iot  directly  related  to  Lagrange  reconstruction. 
Throughout  this  section,  if  is  an  admissible  stencil  for  degree  n,  the  symbol 
denotes  the  convex  hull  of  the  union  of  the  elements  of 


2.2.1  The  case  of  a  smooth  function 

In  [12]  we  have  shown  the  following  result.  Its  proof  follows  easily  from  Ciarlet  &  Raviart’s 
proof  on  Lagrange  and  Hermite  interpolation. 

Theorem  2.1  Let  S  be  an  admissible  (for  degree  n)  stencil  o/IR^,  let  h  and  p  be  respectively 
the  diameter  of  K{S)  and  the  supremum  of  the  diameters  of  the  circles  contained  in  I\{S). 
Let  u  be  a  function  that  admits  everywhere  in  K{S)  an  n  +  derivative  with 

Mn+\  =  5up{ II u(x) II;  X  €  <  +oo. 


If  Pu  is  the  solution  of  problem  V,  then  for  any  integer  rn,  0  <  m  <  n, 
sup{\\D^u{x)  -  D^P^{x)\\-,  X  €  K{S)}  <  CMn+,^ 

^because  we  have  cissumed  an  algebraic  dependency  of  the  control  volumes  in  terms  of  the  points  of  Ad 
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for  some  constant  C  =  C(m,n,<S).  Moreover,  if  S'  is  obtained  from  S  by  an  affine  transfor¬ 
mation,  that  is  there  exists  Xq  €  IR^  and  an  invertible  matrix  A  such  that 

C'l^  E  S'  iff  there  exists  Ck  E  S  such  that  =  A  Ck  +  xq  ), 

then 

C{n,m,S)  =  C{n,m,S'). 

This  result  basically  states  that  if  the  stencil  S  is  not  too  flat,  i.e.  the  ratio  hj p  is  not  too 
large,  then  Pu  will  be  a  good  approximation  of  u.  Let  us  turn  now  to  the  case  of  nonsmooth 
functions. 

2.3  The  case  of  a  nonsmooth  function 

We  only  discuss  the  case  of  piecewise  smooth  functions.  This  class  is  large  enough  for  our 
purpose.  To  do  the  analysis  we  have  to  introduce  the  following  property  which  prevents 
geometrical  degeneration: 

Property  2.2  Let  us  give  e  >  0.  The  admissible  stencil  5^"^  belongs  to  P"  if  and  only  if: 
for  any  vector  U  whose  components  are  zeroes  and  ones  and  both  values  are  represented,  and 
for  P  the  order  polynomial  defined  by 

for  all 

the  following  property  holds:  the  sum  of  the  absolute  value  of  the  n-th  order  coefficients  of  P 
is  greater  or  equal  to  e,  that  is 

det{Ro  ■  ■  •  •  •  •  Rl  ■  ■  ■  ^ 

deti^Ro  •  ■  ■  /2yv(n-l)  '  '  '  '  '  ‘  ^A/(n) 

where  we  have  adopted  the  lexicographic  ordering  for  the  monomials  so  that 

Rk  stands  for  the  column  of  the  determinant  (3)  and  Rt  =  U. 

Then  we  can  prove  [12]  the  following  theorem  that  describes  the  asymptotic  behavior  of 
the  leading  coefficients  of  the  approximation  in  the  mean  of  a  piecewise  smooth  function: 

Theorem  2.3  Let  e  be  a  positive  real  number  and  S  an  admissible  stencil  for  degree  n  such 
that  there  exists  an  affine  transformation  A  as  in  theorem  2.1  for  which  A{S)  E  V'f .  Let 
{xQ,yo)  be  any  point  of  the  set  K{S)  and  u  be  a  real  valued  function  defined  on  a  open 
subset  Cl  of  containing  R{S).  We  assume  that  u  is  in  (7^“’,  p  <  n,  on  Cl  and,  except 
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on  a  locally  (7*  curve,  admits  a  continuous  and  bounded  derivative  with  a  jump  [D^u], 

|[Z)’’u]|  >  My  >  0.  Then,  the  highest  degree  coefficients  of  the  Taylor  expansion  of  satisfy 

\ai,\>C{n,p,e)-^  (5) 

j+j=n 

where  C{n,p,e)  is  a  constant  independent  of  S  and  invariant  by  affine  transformation. 

2.4  Study  of  the  linear  problem  to  solve  for  the  reconstruction 

In  this  section,  we  study  the  numerical  system  one  solves  to  get  from  the  data.  We 
consider  two  kinds  of  expansions  of 

(i)  The  “natural”  expansion:  for  any  point  (a:o,yo)  €  IR^, 

Pu=  Y.  a,,{X  -  xoYiY  -  y^y  (6) 

*+^'<71 

(ii)  An  expansion  using  “barycentric”  coordinates  that  we  now  describe:  let  =  jC] ,  C2,  C3, . . . , 
be  an  admissible  stencil.  Hence,  at  least  one  subset  of  three  elements  of  5^”^  is  an  ad¬ 
missible  stencil  for  n  =  1.  We  may  assume  that  the  set  {C\,C2-,Cz}  is  admissible.  We 
consider  the  three  polynomials  A,  of  degree  1  defined  by: 

<  Ai  >c,=  1  <  f  <  3,  1  <  i  <  3.  (7) 

The  symbol  Sf  is  the  Kronecker  symbol.  Clearly,  we  have 

Ai  +  A2  +  A3  =  1. 

These  polynomials  are  the  barycentric  coordinates  of  the  triangle  constructed  on  the 
gravity  centers  of  C\,  C2,  and  C3.  In  order  to  get  expansion  6,  a  strategy  may  be  to 
look  first  for  the  expansion  of  the  polynomial  P„  in  terms  of  power  of  A2  and  A3: 

P=  Y.  (8) 

i+j<n 

and  then  to  get  the  Taylor  expansion  of  P„  around  the  barycenter  of  C\  from  (8)  (the 
theorems  2.1  and  2.3  give  the  behavior  of  the  leading  coefficients  of  Pu  whatever  the 
point  chosen  in  the  convex  hull  of  5). 

In  order  to  get  the  expansions  (6)  or  (8),  one  has  to  solve  linear  N{n)  x  N{n)  systems: 

B{aQo  •  •  •  flon)^  =  {<u  >c,.  <u  >c.„,  y  (9) 

•  /V(  n) 

vhere  the  matrix  B  is  obtained  by  taking  the  average  of  (A  —  xo)’(V^  —  1/0 V  for  (6)  and  A2A3 
or  (8).  Let  us  now  study  the  properties  of  the.se  linear  systems. 


8 


2.4.1  Case  of  expansion  (6) 

A  very  easy  consequence  of  the  inequality  (5)  is  that: 

Proposition  2.4  Let  us  assume  that  the  conditions  of  theorem  2.3  holds,  and  let  h  be  the 
supremum  of  the  diameters  of  the  spheres  containing  Then  the  condition  number 

of  system  (9)  is  at  least  0(/i“")  for  h  sjnall  enough. 

Proof:  For  the  sake  of  simplicity  we  consider  the  following  norm  on  1R„[X,  K]:  for  P  = 
^i+j<n  -  ^o)'(V''  -  yoY,  ||P|1  =  Ei+j<7i  Oo  1R^^”\  we  consider  the  L'  norm.  Let 

U  be  a  set  of  data  for  the  right  hand  side  of  (9),  and  consider  the  perturbation  6U , 

SIJ  = 

where  e  is  at  the  /  th  position,  /  >  N{n  —  1)  +  I.  All  the  other  entries  of  U  are  zero.  If  one 
considers  the  function  u  defined  on  \jCi  by 


X  e  Ci,  u{x)  =  SUi, 

one  tail  apply  theorem  2.3.  Hence,  the  perturbation  SP  has  a  norm  satisfying 

II«C||  >  E  l«»ul  > 

i+j=n  " 

since  l|^f/||  =  e.  This  complete  the  proof.  □ 

This  fact  is  well  known  for  ID  Lagrange  interpolation  and  has  motivated  the  search 
for  more  efficient  algorithms,  such  as  the  Newton  algorithm.  There  exist  algorithms  that 
generalizes  it  [18,  19].  They  involve  numerous  solutions  of  linear  systems,  so  that  we  have 
preferred  a  more  classical  approach  (see  section  2.5),  for  which  the  coefficients  of  the  linear 
systems  are  obtained  from  the  “barycentric”  roordinate  expansion  (8)  as  it  is  explained  now. 

2.4.2  Case  of  expansion  (8) 

In  the  case  of  expansion  (8),  we  have  the  following  result: 

Proposition  2.5  If  property  2.2  holds  for  some  t  >  0,  then  the  condition  number  of  the 
system  (9)  for  the  expansion  (8)  is  bounded  above  and  below  by  constants  independent  of  h, 
the  supremum  of  the  diameters  of  the  circles  containing  K{S^^^). 

Proof:  The  proof  is  also  based  on  that  of  theorem  2.3.  As  in  proposition  2.4,  the  only  thing 
that  we  have  to  do  is  to  study  the  effect  on  the  Oifs  of  a  perturbation  SU.  We  denote  by  P 
the  polynomial  whose  averages  are  defined  by  6U,  The  proof  can  be  achieved  in  two  stages: 
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(i)  Let  B  any  invertible  matrix.  Consider  the  stencil 

]  +  a:o| 

t  '■  h<j<N{7i) 

for  any  xq.  It  is  clear  from  the  definition  of  the  A,’s  that  Ai{x)  =  A,(i)  if  5  = 
Bx  +  aro-  Hence,  the  sum  S{P)  of  the  absolute  values  of  the  coefficients  of  P  in  the 
liMsis  A2(x)A3(x)  is  the  same  as  that  of  the  development  of  P  in  the  basis  A2(x)A3(x). 
1  iiis  is  a  homogeneity  property. 

(ii)  Since  the  set  of  stencils  defined  by  property  (2.2)  is  compact,  S{P)  is  bounded  below 
and  above,  independently  of  B,  hence  independently  of  h\ 

Ci>F>  C2{t,n)  >  0. 

The  constant  C'2(e,n)  is  larger  than  zero  because  8U  /  0 
This  achieves  the  proof  □ 


2.5  The  explicit  calculation  of  the  reconstruction 

From  the  previous  results,  the  evaluation  of  the  coefficients  a,j  in  (6)  is  done  through  those 
of  (8)  and  hierarchally.  For  the  sake  of  simplicity,  we  assume  that  for  any  p  <  n,  the  set 
of  the  N{p)  first  elements  of  is  an  admissible  set  for  order  p.  This  can  be  achieved  with 
a  suitable  numbering  of  the  elements  of  The  idea  is  that  instead  of  looking  directly 

for  the  coefficients  of  P^^\  to  get  first  those  of  all  of  the  P^'^^'s,  the  reconstruction  over 
for  1  <  ^  <  n.  Then  to  construct  those  of  In  the  ENO  algorithm  described  in  section 
2.6,  this  involves  no  extra  cost  and  simplifies  the  evaluation  of  the  a^j’s.  This  also  has  the 
advantage  of  reducing  the  size  of  the  linear  systems  and  also  of  improving  their  condition 
numbers. 

Assume  that  P^^\  . . . ,  are  known.  We  first  compute  the  coefficients  of  —  P^^\ 


In  equation  (10),  ^  (respectively  02)  stands  for  the  coefficients  {«ij  }i+j<p  (resp.  ). 

The  block  matrices  Ap,  Bp  p+i,  Cp  p+i  and  Dp  p^.i  are  defined  according  to  this  decomposi¬ 
tion.  In  particular,  we  notice  from  the  hypothesis  that  Ap  is  invertible. 
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From  the  conservation  property,  we  get  ^  =  0,  so  that  the  system  (10)  can  be  split: 


—  ~-^p  Bp  0+1^2 


'Cp  Bp  p^\  "h  D 


p  p+i 


Cp  p+i-'^p  Bp  p+i  “h  Dp  p^i 


Since  is  admissible,  Ep  = 

can  get  then  and  finally  the  coefficients  of 
Simple  manipulations  show  that 


(11) 

02 

is  also  invertible,  so  that  one 


/l-i  - 

^P+\  ~ 


A:^B, 


p  p+t  Bp  *  C, 

-E-^a 


p  p+i^p’  +  Ap^ 


1-1 


p  p^l/lp 


Ap  Bp  p+j 


b;^ 


so  that  one  can  quite  easily  go  to  the  next  step.  In  our  case,  since  the  total  degree  of  the 
reconstruction  is  less  or  equal  to  4,  at  most  two  stages  of  this  method  are  needed. 

Finally,  one  must  notice  that  the  condition  number  of  this  method  is  always  better  that 
that  of  the  original  one  because  it  depends  only  on  part  of  the  original  system. 


2.6  The  E.N.O.  reconstruction 

In  [12],  we  have  found  that  only  a  few  stencils  were  indeed  necessary  to  achieve  an  essentially 
non-oscillatory  reconstruction  of  a  piecewise  smooth  function.  This  set  has  to  be  as  isotropic 
as  possible.  Moreover,  the  ENO  reconstruction  was  found  to  achieve  the  expected  order  of 
accuracy  for  smooth  functions,  even  on  very  irregular  meshes.  In  what  follows  Oij  always 
stands  for  any  of  the  coefficients  of  the  reconstruction  P  in  the  natural  basis,  {(X  —  a:o)‘(T  — 
i/o)-’}.  Let  us  describe  our  procedure  up  to  fourth  order: 

(i)  Let  us  start  from  a  given  cell,  Co  assigned  to  a  point  of  A4,  say  [xo^yo). 

(ii)  Consider  all  the  triangles  having  (a:o,  j/o)  as  a  vertex,  and  choose  the  one,  say  T-min-, 
that  minimize 

Z!  l«ol- 

«+j=i 

Here,  is  the  set  of  control  volumes  corresponding  to  the  vertices  of  Tmin-,  (see 
Figure  3-a).  For  a  regular  unstructured  mesh,  there  are  about  six  possible  triangles. 

(iii)  Consider  Tmm-  For  any  of  its  three  edges,  consider  the  three  triangles,  Ti,  Ti,  Ts  as  in 
Figure  3-a.  There  are  three  possible  configurations.  We  choose  the  one  that  minimizes 
the  sum 

Z 

«-Hj=2 
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(iv)  Consider,  as  in  Figure  3-b,  the  configuration  for  third  order.  It  is  obtained  as  fol¬ 
lows:  for  a  stencil  made  of  the  control  volumes  associated  with  the  vertices  of 
{TminiT\,T2,T3}^  oue  may  consider  its  “edges”  made  of  the  external  sides  of 
{ J2,  Tmin},  {Tmin,  Ts}-  Consider  one  of  them,  say  {T2,  T3},  and  the  vertices  a,  fS  and  7. 
Since  the  triangulation  is  conformal,  there  exist  a  triangle  T4  ^  T2  on  the  other  side  of 
Similarly,  Ts  for  [/i, 7].  Then  one  can  construct  Te,  TV,  Tg  and  Tg,  analogously. 
The  stencils  for  fourth  order  reconstruction  are  the  union  of  <5^^^  and  the  control  vol¬ 
umes  associated  with  the  additional  vertices  of  either  {T4,  T5,  Tg,  T7}  or  {T4,  T5,  Tg,  Tg} 
or  {T4, T5, T7, Tg).  For  a  stencil  there  are  at  most  12  stencils  for  fourth  order 

reconstruction. 

The  situations  seem  to  be  become  more  and  more  complicated  as  the  degree  increases. 
Nevertheless,  there  is  a  very  easy  way  to  simplify  it,  so  that  at  each  level  only  three  new 
stencils  for  the  n  -f  order  have  to  be  considered  over  those  of  a  order  one,  just  as 
from  second  order  to  third  order.  Assuming  a  mesh  Ad,  we  want  to  derive  a  -f  order 
reconstruction  method.  The  idea  is  to  work  with  the  control  volumes  defined  for  a  mesh  Ad', 
the  points  and  the  triangulation  of  which  are  constructed  from  those  of  M  by  adding,  for  each 
triangle  of  Ad,  the  points  and  the  triangles  associated  with  the  Pk  Lagrange  interpolation 
[20]. 


3  A  class  of  high  order  numerical  scheme  for  com¬ 
pressible  flow  simulations 


3.1  The  Euler  equations 

Let  us  quickly  recall  elementary  things  about  the  Euler  equation  of  a  calorically  perfect  gas: 

dW  ^  dF{W)  ^  dG{W)  ^ 
dt  Ox  dy 

As  usual,  in  equation  (12),  W  stands  for  the  vector  of  conserved  quantities  and  F  (respec¬ 
tively  G)  is  the  flux  in  the  x  direction  (resp.  y  direction): 


(  ^  \ 

1 

(  pu  \ 

^  p  u  \ 

w  = 

p  u 

p  V 

F{W)  = 

p  +  p 

p  uv 

GiW)  = 

p  uv 
p  v'^  +  p 

\  E  ) 

\  u{E  +  p)  ) 

\  v{E  +  p)  ^ 

(13) 


with  initial  and  boundary  conditions.  In  equation  (13),  p  is  the  density,  u,v  are  the  com¬ 
ponents  of  the  velocity,  E  is  the  total  energy,  and  p  the  pressure,  related  to  the  conserved 
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quantities  by  the  equation  of  state: 


P  =  (7  -  1) 


1 

2 


p{v?'  + 


>) 


(14) 


The  ratio  of  specific  heats,  7,  is  kept  constant. 

It  is  well  known  that  the  system  defined  by  equations  (12),  (13)  and  (14)  is  hyperbolic: 
for  any  vector  n  =  {nx,ny),  the  matrix 


,  OF  dG 


(15) 


is  diagonalizable  and  has  real  eigenvalues  and  eigenvectors.  Let  us  describe  now  the  con¬ 
struction  of  a  order  scheme. 


3.2  Finite  volume  formulation 

We  consider  a  mesh  M.  and  the  control  volumes  as  in  figure  2.  The  semi  discrete  finite 
volume  formulation  of  (12)  is: 

Here  W{t)\s  the  (spatial)  mean  value  of  W{x,t)  at  time  t  over  Ci,n  =  (n^,  Uy)  is  the  outward 
unit  normal  to  dCi,  and  =  nxF  +  UyG.  We  first  describe  the  spatial  approximation  of 
(16),  then  the  temporal  discretization  of  the  resulting  set  of  ordinary  differential  equations. 
Finally,  we  detail  the  boundary  conditions. 

3.2.1  Spatial  discretization 

For  the  sake  of  simplicity  we  define  the  integer  number  p  such  that  either  k  =  2p  or  k  =  2p+ 1. 
The  first  step  is  to  discretize  Ci{t)  up  to  k‘^  order.  First,  we  can  rewrite  area{Ci)Ci{t)  as 

/  F4W{x,t)]dl  =  '£  f  Fn[W{x,t)]dl  (17) 

J dc,  P  Jr , 

where,  as  in  figure  2,  the  set  Fj’s  is  that  of  the  linear  edges  of  C,.  On  each  Fj,,  n  is  constant. 
We  consider  on  any  the  p  Gaussian  points  {Gi}i<i<p  associated  to  the  Gaussian  formula 
of  order  2p  +  1.  The  integral  /p^  Ffi[W{x,t)]dl  is  approximated  by 

/=1 

where  term  is  defined  now.  Set  Cj  is  the  other  control  volume  of  which  F,  is  a  part 

of  the  boundary.  In  Ci  and  Cj  ,  one  computes  the  ENO  reconstructions  at  time  t  of  W, 
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/?,[^(  •  and  Rj[W{  .  up  to  order  k.  The  ENO  reconstruction  of  section  2  is 

applied  to  the  physical  variables,  then  one  deduce  the  conserved  ones.  From  that,  we  set,  in 
equation  (18): 

QkjW  =  .  ,i),k]{G,),R,lW(  .  ,t),k](G,)] .  (19) 

In  equation  (19),  may  be  any  of  the  available  Riemann  solvers.  In  all  the  example 

below,  we  have  chosen  Roe’s  Riemann  solver  with  the  Harten-Hyman  entropy  correction. 


3.2.2  Temporal  discretization 


The  equations  (16),  (17),  (18)  and  (19)  define  a  finite  set  of  ordinary  differential  equations 
that  we  write  as: 

^TF,(t)  =  £.(t)  (20) 

In  (20),  Ci{t)  is  the  discrete  version  of  £,(<).  This  equation  is  discretized  by  the  order 
version  of  the  Runge-Kutta  scheme  of  Shu  [9]: 


r(m)s 


+  /3, /  =  1, 2, ■■■,?,  £<”■'  = 


(21) 


The  order  of  accuracy,  as  well  as  its  TVD  properties,  is  achieved  by  adequate  sets  of  coeffi¬ 
cients  aim,  l^im,  and  p  (see  [9]  for  details). 


3.2.3  Boundary  conditions 

Let  r  be  the  boundary  of  the  computational  domain  and  n  be  the  outward  normal  unit  on 
r.  We  assume  that  F  is  divided  into  two  parts,  F  =  FonLoo,  on  which  different  boundary 
conditions  will  be  used.  Here,  Fo  represents  a  solid  wall,  while  Fqo  represents  the  far-field 
(inflow  or  outflow). 

We  do  not  treat  a  boundary  conditions  by  forcing  the  value  of  a  variable  to  a  prescribed 
boundary  value,  but  consider  instead  the  integral  formulation  (16)  and  apply  boundary 
condition  by  modifying  the  flux  integrals  on  dCi  for  those  cells  with  rC\dCi  ^  0. 

For  example,  for  a  vertex  i  located  on  Fq,  we  do  not  impose  the  slip  condition  U  •  n  =  0 
but  take  this  condition  into  account  in  the  evaluation  of  the  convective  flux 

/  Fux  -b  Gn 


0 

fpo  n  P^y 
0 
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For  a  vertex  located  on  Too,  we  again  use  an  approximate  Riemann  solver.  We  define  a 
far  field  state  Woo-,  =  fraof)9C,  ^  agreement  with  what  has  been  done  in  the 

interior  of  the  computational  domain, 


/  Fn^  +  Gnj,  =  4>(W,-,  VFoo,n.).  (22) 

2roon9‘'’> 

In  equation  (22),  is  a  numerical  flux  function.  For  simplicity,  we  have  chosen  a  modified 
Steger- Warming  flux  instead  of  the  Roe  flux. 


^W„  Woo,n)  =  A+W,  + 

The  matrices  At  and  Az  are  the  positive  and  negative  parts  of  the  matrix  An  defined  in 
(15)  and  evaluated  for  W  =  W^. 

In  all  the  examples  we  have  treated  below,  the  boundaries  were  either  fully  subsonic 
or  fully  supersonic,  so  that  the  procedure  was  really  simple,  contrary  to  what  would  have 
appeared  in  a  mixed  type  boundary  condition. 

Finally,  we  have  reduced  the  order  of  accuracy  of  the  reconstruction  for  cells  that  are  too 
near  to  the  boundary.  For  them,  a  proper  calculation  of  the  ENO  stencil  may  be  impossible 
because  the  set  of  possible  stencils  is  biased  in  one  direction  due  to  the  boundary.  For  the 
third  order  scheme,  these  cells  are  those  related  to  a  mesh  point  that  belongs  to  a  triangle 
having  at  least  one  point  on  the  boundary.  For  the  fourth  order  scheme,  they  are  those 
belonging  to  a  triangle  which  shares  a  vertex  with  a  cell  for  which  a  reduction  of  order  must 
be  done  for  third  order. 

3.3  Positivity  of  the  density  and  the  pressure 

As  pointed  out  by  Harten  et.  al.  [22],  in  some  situation  and  for  a  very  few  cells,  the  ENO 
reconstruction  of  the  density  and  pressure  may  lead  to  negative  values.  For  these  cells,  and 
these  cells  only,  following  [22],  we  reduce  the  order  of  accuracy  with  the  following  inductive 
method  {w  is  either  the  density  or  the  pressure,  w,  is  its  average  on  C,).  Consider,  in  C,,  the 
reconstruction 

HKn|(A:,r)  =  y;  a„(,x  -  x,nY  -  y,)\ 

1=0  p+q=t 

If  Yl?=2^p+q=t  ~  ^gY{y  ~  ygY\  ^  at  a  Gaussian  point  (x,y),  then  the  recon¬ 

struction,  for  that  point,  is  set  to  R[w,n  —  l](a:,i/).  Then,  we  repeat  the  test  if  necessary. 
Usually,  the  parameter  a  is  set  to  0.95. 


In  all  the  simulations  we  have  done,  the  tests  were  positive  for  a  very  small  set  of  cells 
and  a  zeroth  order  reconstruction  was  never  used.  They  were  never  positive  for  the  second 
order  scheme.  The  number  of  points  causing  difficulty  is  problem  dependent.  For  example, 
only  three  points  caused  occasional  problems  for  the  facing  step  problem  with  a  5000  node 
mesh.  These  points  were  all  located  in  the  front  shock. 


4  Numerical  tests 

All  the  example  we  propose  now  have  been  computed  with  the  second  and  third  ENO 
schemes.  The  ratio  of  specific  heats,  7  is  always  set  to  1.4. 

4.1  A  linear  advection  problem 

In  order  to  test  the  precision  of  these  scheme,  we  have  computed  the  advection  of  sine  wave 
on  a  sequence  of  totally  unstructured  meshes  with  an  increasing  number  of  points.  The  con¬ 
vection  velocity  was  parallel  to  the  x  axis  but  since  the  meshes  was  totally  unstructured,  this 
is  not  a  privileged  direction.  Figure  4  shows,  in  the  abscissa,  the  logarithm  of  the  maximum 
radius  of  the  circumscribed  circles  of  the  triangles  of  the  meshes,  and  in  the  ordinate,  the 
logarithm  of  the  maximum  absolute  value  of  the  difference  between  the  computed  values. 
The  exact  solution  is  also  indicated.  The  slopes  —2  and  —3  are  indicated,  so  that  one  can 
see  that  the  expected  order  of  accuracy  is  indeed  achieved.  Figure  5  shows,  for  the  medium 
mesh,  a  cross  section  of  the  computed  solution  for  the  second  and  third  order  scheme.  The 
exact  cross  section  is  also  indicated.  One  can  see  that  the  main  differences  lies  at  the  extrema 
of  the  sine,  as  expected. 


4.2  A  Shock  tube  problem 

We  have  set  up  a  two  dimensional  shock  tube  problem  in  the  square  [0, 1]  x  [0, 1].  Its 
boundary  are  solid.  The  initial  conditions  are: 


f  ^=1.0 

forx  <  0.5  and  \y  —  0.5|  <  0.25,  <  u  =  v  =  0.0 

[  p=\.0 


else 


< 


p  =  0.125 

u  —  V  =  0.0 
p  =  0.1 


The  mesh  is  completely  unstructured  with  2127  nodes  and  4088  triangles.  The  velocity  field 
obtained  by  the  third  order  scheme  at  time  t  =  0.9  is  displayed  in  Figure  6.  The  differences 


16 


between  both  results  are  more  clearly  visible  in  the  near  stagnation  zone.  In  order  to  better 
represent  that  area,  we  have  removed  from  the  velocity  field  all  the  points  for  which  the  sum 
of  the  absolute  values  of  the  two  components  is  larger  than  0.15.  The  result  is  shown  in 
Figures  7  (second  order)  and  8  (third  order).  One  can  clearly  observe  that  the  number  of 
small  structures  of  the  flow  is  much  more  important  in  Fig.  8  than  in  Fig.  7.  The  shocks  in 
the  upper  and  lower  part  of  the  pictures  are  also  resolved  differently.  Their  location  is  also 
different,  though  this  can  be  seen  only  by  superimposing  the  pictures. 

One  should  also  mention  that  this  test  is  not  particularly  easy  for  our  method.  After  a 
short  time,  the  shock  reflects  from  the  wall.  The  reflected  shock  interacts  with  the  others 
structures  of  the  flow,  leading  to  interactions  between  the  various  kind  of  discontinuities 
and  with  the  smooth  parts  of  the  flow.  The  multiple  regions,  as  shown  on  our  figures,  with 
different  kind  of  discontinuities  (contact  and  shock)  are  resolved  by  our  method  without  any 
special  tricks. 

4.3  A  Mach  3  wind  tunnel  with  a  step 

We  have  run  this  test  case,  documented  in  [23],  for  the  second  order  and  third  order  ENO 
schemes  on  a  5140  node,  9958  triangle  mesh.  This  discretization  corresponds  to  the  medium 
mesh  used  in  [23].  A  portion  of  it  is  displayed  in  Figure  9.  It  is  totally  unstructured.  The 
conditions  of  the  problem  are  the  following:  a  uniform  Mach  3  flow  is  set  in  a  channel.  At 
the  initial  time,  a  step  of  relative  height  0.2  is  installed  in  the  channel.  The  channel  length 
is  3  and  the  step  is  located  at  0.6.  This  situation  creates  a  shock  that  reflects  on  the  upper 
part  of  the  channel  then  evolves  to  a  lambda  shock  as  time  increases.  It  interacts  with  the 
upper  part  of  the  step.  A  weak  shock  is  also  created  by  the  expansion  wave  at  the  corner. 
This  shock  interacts  with  the  reflected  one  creating  a  slip  line.  The  location  of  this  slip  line 
is  highly  dependent  on  the  boundary  conditions  set  at  the  corner. 

Here,  no  special  treatment  is  done,  contrary  to  what  was  advocated  in  [23],  so  that 
the  quality  of  the  second  reflected  shock  is  poor.  We  only  want  to  verify  the  effect  of  the 
increasing  order  of  accuracy  on  the  solution,  so  that  we  will  only  look  at  the  first  reflected 
shock.  The  solutions  of  Figure  10  (second  order  ENO)  and  Figure  11  (third  order  ENO)  are 
shown.  A  clear  improvement  on  the  thickness  of  that  reflected  shock  can  be  seen  from  the 
horizontal  cross  section  of  the  density  at  y  =  0.5,  Figure  12.  The  slip  line  coming  from  the 
lambda  shock  is  also  more  visible  in  Figure  11  than  in  Figure  10  as  well  as  the  weak  shock 
near  the  corner. 
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4.4  Reflection  of  a  shock  on  a  wedge 

This  problem  is  also  well  documented  in  the  literature.  In  order  to  achieve  a  correct  solution, 
one  has  either  to  use  very  fine  meshes  or  adapted  meshes  (see  [5]  for  example).  We  have 
chosen  a  case  where  the  planar  shock  enters  from  the  left  in  a  quiescent  fluid.  Its  Mach 
number  is  Ms  =  5.5  and  is  defined  for  the  flow  values  in  the  quiescent  fluid  where  the 
density  is  set  to  1.4  and  the  pressure  to  1.  One  expects  a  double  Mach  reflection. 

The  mesh  has  only  8569  points  and  16806  triangles.  A  part  of  it  is  shown  in  Figure  13. 
The  density  contours  of  the  two  calculations  are  displayed  in  Figures  14  (second  order)  and 
Figure  15  (third  order).  A  very  clear  improvement  of  the  slip  line  coming  from  the  Mach 
stem  can  be  observed.  The  second  triple  point  can  also  been  observed  in  Figure  15,  though 
it  is  of  poor  quality  because  of  the  insufficient  resolution  of  the  present  mesh.  It  is,  however, 
totally  obscured  in  Figure  14.  Generally  speaking,  all  the  discontinuities  are  better  resolved 
by  the  third  order  scheme. 

5  Conclusions 

A  third  order  ENO  scheme  has  been  derived  for  triangular  unstructured  meshes;  this  demon¬ 
strates  the  possibility  of  deriving  ENO  schemes  for  unstructured  meshes.  We  indicate  how 
to  build  higher  order  ENO  schemes  and  give  some  comments  on  the  numerical  stability  of 
the  reconstruction  step. 

Our  new  scheme  has  been  tested  on  a  set  of  well  known  test  cases  and  compared  to 
a  second  order  scheme.  In  all  cases,  the  results  are  clearly  improved.  Our  results  also 
demonstrate  the  scheme’s  robustness.  The  cost  of  the  scheme  is  four  times  that  of  the  second 
order  scheme  (on  a  Cray  YMP)  even  though  the  new  code  is  far  from  being  optimized.  In 
particular,  no  optimization  has  been  done  in  the  ENO  reconstruction  procedure,  the  most 
expensive  routine,  so  that  the  factor  of  four  is  clearly  a  poor  upper  bound  on  the  cost  ratio. 

In  the  near  future,  we  will  derive  the  fourth  order  version  of  this  class  of  schemes.  The 
two  schemes  will  then  be  coupled  with  a  dynamic  adaptation  procedure  [3]  to  improve  their 
efficiency. 
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Figure  1:  Covering  of  the  rectangle  [xq,  Xj]  x  [yo,  yi]  by  trieingular  control  volumes 
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Edge  of  the  control  volume 
Figure  2:  Element  of  the  dual  mesh 
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Figure  6:  Shock  tube  :  velocity  field  at  time  t  =  0.9 
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Figure  7:  Zoom  of  the  velocity  field  in  [0.5,1]  x  [0.25,0.75].  Second  order  solution 


Figure  8:  Zoom  of  the  velocity  field  in  [0.5, 1]  x  [0.25, 0.75].  Third  order  solution 
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Figure  10;  Density  contours  for  the  second  order  ENO  solution,  t=4,  min=0.329,  max=4.64. 
Density  contours  from  0.287  to  4.584,  A/j  =  0.14 


Figure  11:  Density  contours  for  the  third  order  ENO  solution,  t=4,  min=0.287,  max=4.584, 
A/»  =  0.14 
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Figure  13;  Portion  of  the  mesh  used  for  the  Mach  reflection  problem 


Figure  14;  Refection  of  a  planar  shock  by  a  wedge  :  density  contours,  second  order  solution. 
Min=1.4,  Max=17.3.  Contour  from  1.4  to  19.088,  Ap  =  0.36 
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